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An  expansion  wave  is  produced  when  an  incident  shock  wave  interacts  with  a  surface  separating  a  fluid  from 
a  vacuum.  Such  an  interaction  starts  the  feedout  process  that  transfers  perturbations  from  the  rippled  inner 
(rear)  to  the  outer  (front)  surface  of  a  target  in  inertial  confinement  fusion.  Being  essentially  a  standing  sonic 
wave  superimposed  on  a  centered  expansion  wave,  a  rippled  expansion  wave  in  an  ideal  gas,  like  a  rippled 
shock  wave,  typically  produces  decaying  oscillations  of  all  fluid  variables.  Its  behavior,  however,  is  different  at 
large  and  small  values  of  the  adiabatic  exponent  y.  At  y>3,  the  mass  modulation  amplitude  Sm  in  a  rippled 
expansion  wave  exhibits  a  power-law  growth  with  time  ocf£  where  /3=(y- 3)/(y- 1).  This  is  the  only  example 
of  a  hydrodynamic  instability  whose  law  of  growth,  dependent  on  the  equation  of  state,  is  expressed  in  a  closed 
analytical  form.  The  growth  is  shown  to  be  driven  by  a  physical  mechanism  similar  to  that  of  a  classical 
Richtmyer-Meshkov  instability.  In  the  opposite  extreme  y-l<Sl,  Sm  exhibits  oscillatory  growth,  approxi¬ 
mately  linear  with  time,  until  it  reaches  its  peak  value  ~(y-  1)_1/2,  and  then  starts  to  decrease.  The  mechanism 
driving  the  growth  is  the  same  as  that  of  Vishniac’s  instability  of  a  blast  wave  in  a  gas  with  low  y.  Exact 
analytical  expressions  for  the  growth  rates  are  derived  for  both  cases  and  favorably  compared  to  hydrodynamic 
simulation  results. 
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I.  INTRODUCTION 

Linear  stability  analysis  of  rippled  centered  rarefaction 
waves  has  been  performed  in  recent  years  in  connection  with 
the  studies  of  the  classical  Richtmyer-Meshkov  (RM)  insta¬ 
bility,  reflected  rarefaction  case  [1-4].  Rippled  rarefaction 
waves  were  found  to  be  unstable.  Instability  of  a  reflected 
rippled  rarefaction  wave  produced  in  a  shock-interface  inter¬ 
action  is  manifested  in  the  linear  growth  of  the  ripple  ampli¬ 
tude  at  its  trailing  edge  [2-4].  The  ripple  amplitude  of  its 
leading  edge  remains  constant  in  time.  The  unstable  trailing 
edge  of  a  rippled  rarefaction  wave  is  a  weak  discontinuity: 
The  density  gradient  is  discontinuous  there.  This  is  why  this 
instability  of  a  rippled  rarefaction  wave  in  some  sense  also 
turns  out  to  be  weak:  The  linear  growth  of  its  trailing  edge 
displacement  does  not  translate  into  a  similar  growth  of,  say, 
observable  mass  modulation  amplitude  Sm  =  SJ  pdx.  As  first 
shown  in  [3],  for  ideal  gases  with  moderate  values  of  y,  like 
|  or  Sm  in  a  rippled  rarefaction  wave  exhibits  decaying 
oscillations. 

An  expansion  wave  is  a  rarefaction  wave  in  which  a  gas 
rarefies  to  zero  pressure,  expanding  into  vacuum.  Such  an 
expansion  wave  is  produced  in  the  limiting  case  of  a  heavy- 
to-light  shock-interface  interaction,  when  the  density  of  the 
light  fluid  tends  to  zero.  This  case  is  of  interest  for  inertial 
confinement  fusion  (ICF)  studies  because  of  its  relevance  to 
the  so-called  feedout  process  that  transfers  mass  perturba¬ 
tions  from  the  inner  (rear)  to  the  outer  (front)  surface  of  a 
laser-or  hohlraum-driven  target,  producing  a  seed  for 
Rayleigh-Taylor  (RT)  growth  [5-10],  Since  the  small  RT 
seeds  contributed  by  the  outer  and  inner  surface  roughness 
and  the  drive  nonuniformity  are  statistically  independent. 


each  of  these  contributions  could  be  studied  separately,  one 
Fourier-Legendre  mode  at  a  time.  Focusing  on  the  contribu¬ 
tion  of  the  inner  surface  roughness,  consider  a  smooth  shock 
wave  that  arrives  from  the  uniformly  irradiated  outer  surface 
of  the  target  and  breaks  out  at  its  rippled  inner  surface.  Then 
a  rippled  reflected  rarefaction  wave  starts  propagating  back 
to  the  outer  surface,  bringing  mass  perturbations  there.  When 
a  spherical  target  is  empty,  or  filled  with  a  gas  whose  density 
is  much  less  than  the  density  of  its  shell,  or  if  a  planar  target 
in  vacuum  with  ripples  on  its  rear  surface  is  used  to  study 
feedout,  as  in  the  experiments  [7,8,10],  then  the  wave  re¬ 
flected  from  the  inner  (rear)  surface  of  the  target  is  a  rippled 
expansion  wave.  The  limiting  transition  from  a  heavy-to- 
light  case  of  RM  instability  to  the  feedout  situation  when 
there  is  no  light  fluid  behind  the  rippled  surface  reached  by 
the  shock  wave  is  discussed  in  detail  in  [9]. 

As  the  pressure  behind  the  trailing  edge  of  a  rarefaction 
wave  tends  to  zero,  the  growth  rate  of  its  instability  vanishes 
[3,4].  Indeed,  the  fluid  particles  initially  located  at  the  free 
rear  surface  of  the  target,  after  the  shock  interaction  with  this 
surface,  form  the  expansion  front.  In  planar  geometry,  nor¬ 
mal  velocity  of  the  expansion  front  is  the  same  for  its  peaks 
and  valleys,  2as/(y-l),  where  as  is  the  speed  of  sound  be¬ 
hind  the  incident  shock  wave.  Therefore,  the  displacement 
amplitudes  of  both  the  leading  and  the  trailing  edges  of  a 
planar  expansion  wave  are  constant  in  time,  and  in  this  sense 
a  rippled  expansion  wave  is  stable  for  any  y. 

Surprisingly,  it  was  found  in  [9]  that  for  y>  3  the  mass 
modulation  amplitude  Sm  in  a  rippled  expansion  wave  ap¬ 
pears  to  exhibit  an  unlimited  power-law  growth.  In  the  op¬ 
posite  limiting  case,  for  y  close  to  unity,  Sm  is  an  oscillating 
function  of  time,  whose  amplitude  is  also  found  to  grow, 
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linearly  with  time,  until  it  reaches  some  peak  value.  This  new 
instability,  developing  inside  a  rippled  expansion  wave 
bounded  by  stable  leading  and  trailing  edges,  is  thus  unre¬ 
lated  to  the  instability  of  rippled  rarefaction  waves  discussed 
previously  [2-4],  These  findings  could  not  be  fully  under¬ 
stood  at  the  time  because  the  physical  mechanisms  driving 
this  instability  remained  unknown,  and  no  formulas  were 
available  for  its  growth  rate.  Where  does  this  new  instability 
come  from? 

We  report  here  the  solution  of  this  problem.  We  have  ex¬ 
plained  the  physical  mechanisms  driving  the  instability,  de¬ 
rived  explicit  asymptotic  expressions  for  the  growth  rates 
(Sec.  II),  and  favorably  compared  our  theoretical  predictions 
to  2D  numerical  simulation  results  (Sec.  III).  In  Sec.  IV  we 
conclude  with  a  discussion. 


H.  THEORY 

Let  a  planar  slab  of  thickness  L  be  filled  with  a  cold 
uniform  ideal  gas,  with  a  single-mode  ripple  at  the  rear  sur¬ 
face  separating  the  gas  from  vacuum: 


p(x,y,t  — >  -°o)  = 


0, 


x  <  -  L\ 
p0,  -L^x 
0, 


:  Sx() 


:0  cos  ky; 
x  >  Sx o  cos  ky , 


(1) 


where  k=2'irl'K  is  the  ripple  wave  number.  We  assume  the 
ripples  to  be  small  (A:<5x0<81)  and  the  slab  to  be  thick  (kL 
S>  1)  compared  to  the  ripple  wavelength  X..  The  former  con¬ 
dition  allows  us  to  use  the  linearized  small-amplitude  theory. 
Due  to  the  latter  condition,  we  can  follow  the  evolution  of  an 
isolated  rippled  expansion  wave  launched  at  the  rear  surface 
of  the  slab  for  an  arbitrarily  long  time  before  its  leading  edge 
breaks  out  at  the  front  surface.  Mass  modulation  amplitude 
defined  as 

r+<»  r+o° 

Sm(t)  =  S I  pdx  =  |  Spdx 

J  —OO  J  —CO 

1  f+0° 

=  -  [p(x,y  =  0 ,t)~  p(x,y  =  \/2 ,t)]dx  (2) 

^  J  -00 


is  the  integral  measure  of  mass  perturbation  in  the  slab, 


which  is  directly  observable  in  face-on  x-ray  radiographic 
images  in  ICF-related  experiments  [7,10].  Substituting  (1) 
into  (2),  we  find  the  initial  value  of  mass  modulation  ampli¬ 
tude:  Sm0=p0Sx0 ,  the  initial  density  times  the  initial  ripple 
amplitude. 

Let  a  constant  pressure  ps  be  applied  at  t=tx  to  the  front 
surface  of  the  slab.  It  launches  a  strong  planar  shock  wave 
whose  front  velocity,  mass  velocity,  and  the  postshock  speed 
of  sound  are,  respectively. 


D  = 


(y+i  )psYn 

2p0  / 


2 

U= - D, 

y  +  1 


[2y(y-D]1/2 
y  +  1 


-D. 


(3) 


Taking  tx  =  -LID ,  we  make  the  shock  wave  arrive  at  the  rear 
surface  of  the  slab  at  t= 0.  Obviously,  at  f <  0  the  mass 
modulation  amplitude  is  constant,  Sm=Sm0,  because  only 
the  initial  ripples  at  the  rear  surface  of  the  slab  contribute  to 
the  integral  (2).  At  t= 0,  the  planar  shock  wave  breaks  out  at 
the  rippled  rear  surface  of  the  slab,  instantly  producing  a 
rippled  expansion  wave  that  propagates  from  the  rear  surface 
into  the  shocked  gas.  The  reflected  expansion  wave  is 
bounded  by  its  rippled  leading  edge  and  the  expansion  front, 
whose  ripple  amplitudes  are,  respectively, 

8x,=  — '—{y-  1  +  [2y(y-  l)]1/2}&0, 
y+  1 


Sx. 


y+  1 


y-  1  -2 


2y 
y- 1 


1/2 


Sxn 


(4) 


they  remain  constant  in  time,  see  [9].  Starting  from  f=0,  the 
mass  modulation  amplitude  begins  to  change. 

In  what  follows,  we  use  normalized  time 


r  =  kast  = 


27 TXf 


(5) 


where  xt  is  the  distance  traveled  by  the  leading  edge  of  the 
expansion  wave  in  the  shocked  gas.  The  mass  modulation 
amplitude  Sm  will  be  expressed  below  in  units  of  Sm0,  so 
that  Sm{r=  0)=  1.  An  exact  expression  for  Sm( t)  as  a  Taylor 
series  in  t  had  been  obtained  in  [9].  We  present  it  here  in  a 
slightly  modified  form: 


Sm(r)  =  2(3  -  y)2  77 

7-0  J 


"  2(y-  1)  3’ 

{2y'[y-  1  +  V2y(y-  1)]  +  y+  l}r 

j  + 

3  -  y  2 

[  (y—  l);'+  l]{4[(y-  l)y  +  l]2  -  (y—  l)2}r 

y+ 1  1 

j  + 

3  -  y  2 

(6) 


The  infinite  series  (6)  converges  for  any  finite  t,  at  any 
y>  1,  thus  generating  an  entire  analytic  function  of  complex 
variable  t.  Retaining  a  sufficient  number  of  terms  in  the  se¬ 
ries,  we  can  follow  the  evolution  of  mass  modulation  ampli¬ 


tude  for  any  finite  time  rS*  0.  Some  examples  are  shown  in 
Fig.  1. 

Observe  a  somewhat  different  behavior  for  moderate  (a), 
large  (b),  and  small  (c)  values  of  y.  For  the  moderate  y  case 


046306-2 


INSTABILITY  OF  A  PLANAR  EXPANSION  WAVE 


PHYSICAL  REVIEW  E  72,  046306  (2005) 


seen  to  grow  for  a  longer  time  than  in  Fig.  1(a),  and  reach 
higher  peak  values,  exceeding  dm0  by  an  order  of  magnitude 
or  more,  before  the  oscillations  start  to  decay. 

Of  course,  the  time  histories  shown  in  Fig.  1  do  not  con¬ 
stitute  a  proof  of  either  the  decay  of  oscillations  of  dm  in 
Figs.  1(a)  and  1(c)  or  of  the  unlimited  growth  of  dm  in  Fig. 
1(b).  At  some  later  time,  not  shown  in  the  figures,  the  decay 
might  give  way  to  growth,  and  vice  versa.  To  make  an  un¬ 
ambiguous  conclusion,  one  needs  to  determine  the  exact 
asymptotic  behavior  of  the  analytic  function  generated  by  the 
series  (6)  in  the  limits  t— >  °°,  finite  y-  1 ,  and  y—  1  — » 0,  finite 
t,  which  is  done  below  in  this  section.  Once  the  asymptotics 
of  dm  are  established,  we  need  to  understand  them  qualita¬ 
tively:  If  there  is  an  instability,  we  have  to  explain  what 
physical  mechanism  drives  it. 


A.  Power-law  growth  and  oscillations  at  y=(3n-l)/(n+l)  and 

7-0 


At  certain  values  of  y,  the  asymptotic  formulas  for  dm  (  t) 
in  the  limit  t—>° °  could  be  derived  directly  from  Eq.  (6).  To 
do  this,  we  observe  that  in  some  particular  cases  the  right- 
hand  side  of  Eq.  (6)  reduces  to  a  Taylor  expansion  of  a  gen¬ 
eralized  hypergeometric  function.  For  example,  substituting 
y=  |  into  (6),  we  reduce  it  to 


dm  (  t)  =  3E4I 


'  3  3  5  3  1  f\ 

,1’2’2’3’2’4’4’_  8  / 

5  +  \)tl  /  5  5  7  7  5 

- —  ,F4  2,— ,  —  ;4,  — ,  — ,— 

10  3  4\  2  2  4  2  4 


(7) 


FIG.  1.  Time  histories  of  areal  mass  variation  in  a  rippled  ex¬ 
pansion  wave  for  moderate  y=  |  (a),  large  y=3  to  7  (b),  and  small 
y=1.05  (c). 

of  Fig.  1(a),  dm  rapidly  changes  phase,  as  experimentally 
observed  in  [10],  and  then  continues  to  oscillate.  Amplitudes 
of  the  first  several  peaks  of  dm  exceed  its  initial  value  dm0 
by  a  factor  of  3-4;  later,  the  oscillations  decay.  For  a  large  y 
[Fig.  1(b)],  we  notice  a  monotonic  growth  of  dm  superim¬ 
posed  on  the  oscillations:  at  y=3,  dm  does  not  seem  to  decay 
with  time,  whereas  at  y>  3  the  growth  appears  to  be  unlim¬ 
ited.  For  a  small  7  [Fig.  1(c)],  the  oscillation  amplitude  is 


Such  reduction  turns  out  to  be  possible  for  an  infinite 
number  of  values  of  7,  which  satisfy  the  requirement  that  the 
value  of  the  parameter  n  introduced  in  [1,4]: 


3-7 

is  an  integer.  Positive  values  of  1  generate  7=  (3 n 
-  1 )/(«+  1)  between  1  and  3:  1,  v,  2,  y , . . .  .  Negative  values 

13 

of  72^-2  generate  y’s  between  7  and  3:  7,  5,  y,  4,...  .  For 
positive  72  S5  2  the  analytic  function  generated  by  the  series 
(6)  is 


dm  (r)  = 


n+l1  n+2 


1 


72  +  1 


272  -  1  272  -  1  372  -  1  1 


272  -1 


n  -  1  ’  2(72  -  1)  ’  2(72  -  1)  ’  ’2(n-l)  72  —  1  2(t2  —  1)  ’  2t2  ’  2n 

n  372  -  1  2 


1/72-1 
4 


T2 


(72  +  1)(v(3t2  -  1 ) (72  -  1)  +  72  -  IJt2 


2(272  -  1)(3t2-  1) 


n+lr  n+2 


72  -  1  ’  2(n  -  1)  ’  2(72  —  1) 


272  -  1  372 

+  1,  - -  +  1 ; 


■2  572  -  3  1 


2(72  -  1)  72  -  1  ’  2(77  —  1 )  ’  2/2 


+  1, 


2«-  1 


277 


+  1;-- 
4 


1  /  72  —  1 


n-1 


T2  . 


(9) 


For  negative  72  ^  -2 
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Smir) 


/  n  + 1  1  2/7  +  1  1  2/7-13/7-1  1  2«+l  1//7-IW  \ 

"  3_”\«  -  1  ’  2(«  -  1) ’  2h’  ’  2/z  ’  2(/?  —  1)  ’  n  —  1  ’  2{n  —  1) ’  2(/7  —  1 )  ’  ’2(«  —  1)  ’  4 \  n  )  ) 

( n  +  l)(V(3/7  -  l)(/7  -  1)  -  n  +  1)7^  /  77  3/7  -1  1  2/7+1  2/7-13/7-25/7-3 

2(2n  -  1)(3«  —  1)  2_"  3_"\/7  -  1  ’  2(77  -  1) ’  2//  ’  ’  2/?  ’  2(77  —  1) ’  n  —  1  ’2(77- 1)’ 


1  2/7+1  l/77-l 

2(77  -  1)  +  ’  ’  2(/7  -  1)  +  ~  4 


T2  ■ 


(10) 


If  negative  n  is  an  even  number,  then  one  upper  index  [(// 
+  l)/2(/7-l)  in  the  first  hypergeometric  function]  coincides 
with  one  of  the  lower  indices,  reducing  the  order  of  the  hy¬ 
pergeometric  function  to  l_nF2_n. 

The  advantage  of  the  explicit  expressions  (7),  (9),  and 
(10)  is  that  we  can  readily  apply  the  known  formulas  for  the 
asymptotic  expansions  of  the  generalized  hypergeometric 
functions  at  large  negative  argument  [11]: 


pFp+Mi,  ...,ap\bu  ...  ,bp+l;-z2) 


!£>,) 

n:=1  T{ak) 


z2x 

—  cos(2  z  +  ttx) 


+ 


nt^rfa-Q/) 

n;;r  {bj-ad 


(id 


at  z— +00,  where  X=2^k=\ak~ ^j=i^j+V>  and  ab  *  =  1 ,  •••  ,P, 
bj,  j=  1 , ...  ,p+ 1  are  supposed  to  satisfy  the  condition  that 
none  of  ah  bj,  bj-ah  ak-a ,■  is  zero  or  a  negative  integer  for 
any  i,k=\, ...  ,p,  j=  1  ,...,p+ 1,  k=/=i.  Inspecting  Eqs.  (9) 
and  (10),  we  find  only  one  such  case,  when  cij-ak  is  zero:  in 
(9),  if  /7  is  even,  a2=an/2+x  ■  Modification  of  (11)  for  this  case 
is  done  via  a  limiting  transition  a„/2+i— and  yields  two 
terms  in  this  sum,  one  proportional  to  z~2“2  and  another  pro¬ 
portional  to  z~2“2  In  z,  which  is  to  be  expected  in  the  case  of 
degeneracy.  We  will  not  need  the  explicit  form  of  these  terms 
because  they  are  never  dominant  in  the  modified  sum  (11). 
Indeed,  a2=cix  +  \l'2> a.\.  Since 


1  3-7 

/7-1  2(y  —  1) 


(12) 


is  the  lowest  of  all  dj  in  both  cases  (9)  and  (10),  we  note  that 
for  the  hypergeometric  functions  in  (9)  and  (10),  z~2ai  always 
is  the  leading  power-law  term  in  (11).  Substituting  the  values 
of  indices  from  (9)  and  (10)  into  (11),  we  find  that  the  main 
contribution  to  the  oscillating  term  [from  the  second  terms  in 
the  right-hand  sides  of  Eqs.  (9)  and  (10)]  is  determined  by 
X=~2  for  any  n.  Using  the  Gauss’  multiplication  formulas 
([12],  Eq.  6.1.20)  to  simplify  the  products  of  gamma  func¬ 
tions,  and  expressing  the  integer  n  via  the  adiabatic  exponent 
y  in  the  final  results  with  the  aid  of  (8),  we  reduce  the  lead¬ 
ing  terms  of  the  asymptotic  formula  for  Sm  at  r— >°°  to  the 
following  sum: 


Sm  (  t)  s  Smosc(  t)  +  5msec(r) .  (13) 


The  dominant  oscillatory  contribution  coming  from  the  lead¬ 
ing  cosine  term  in  (11)  for  any  y  is  given  by 


Smosc{T)  = 


4 \  2[ y -  1  +  \'2 y(y-  l)]sin  kli 


(y-  l)~flT 


where  17  = 


2(y-l) 
.  y- 1 


(y-i)/(3-y) 


(14) 


in  agreement  with  [9]. 

The  dominant  monotonic  (secular)  contribution  coming 
from  the  leading  power-law  term  in  (11)  for  y<3  describes 
power-law  decay 


<?771sec(7-)  =  - 


2  sinl  - I 

— - - y  J  [2(y-  l)2  +  (y~  3)  \j2y(y-  1)] 

(y- 1) 


y+  1  2y^/2^‘) 
3  -yT  ) 


(15) 


whereas  for  y>3  it  describes  power-law  growth 


Smsec(T)  =  - 


2 

(y-D4 


[2(y-  l)2  +  (y-  3)\;2y(y-  1)J 


y+  1 

y-3 


(r-m(y-i) 


(16) 


Of  course,  (15)  represents  a  branch  of  the  same  analytic 
function  of  r  as  (16).  As  y-3  changes  sign,  the  noninteger 
power  of  negative  (y+ l)r2/(y-3)  in  (16)  becomes  a 
complex-valued  function.  The  right-hand  side  of  (15)  gives 
the  real  part  of  this  complex  function  for  one  of  its  branches. 
The  other  branches  would  produce  phase  shift  factors 
(-1)*  sin[(2i+  l)7r/(y- 1)],  where  s  is  an  arbitrary  integer 
number;  the  branch  (15)  corresponds  to  s= 0,  the  phase  shift 
factor  being  sin[Tr/(y-  1)]. 

The  exact  asymptotic  expressions  ( 1 4)— ( 1 6)  have  been  de¬ 
rived  above  only  for  a  discrete  set  of  y’s  given  by  (8).  Since 
this  is  an  infinite  set  with  an  accumulation  point,  one  can 
reasonably  expect  these  formulas  to  be  valid  for  an  arbitrary 
y>l,  if  the  asymptotic  expressions  sought  for  are  analytic 
functions  of  the  parameter  y  and  pertain  to  the  same  sheet  of 
the  Riemann  surface.  A  rigorous  derivation  of  (16)  in  the 
next  section  confirms  that  this  is  indeed  the  case. 
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FIG.  2.  Exact  (solid)  and  asymptotic  (dotted  lines)  time  histories 
of  areal  mass  variation  in  a  rippled  expansion  wave  for  y=  |  (a), 
y=3  to  7  (b),  and  y=  1.003  (c). 

Asymptotic  formulas  ( 14)— (16)  describe  oscillatory  decay 
and  power-law  secular  growth  of  8m  at  large  t  for  y<  3  and 
y>  3,  respectively.  For  y=3  mass  modulation  amplitude  at 
large  r  tends  to  a  constant  value,  exactly  equal  to  its  initial 
value  but  with  inverted  phase.  Taking  the  limit  y— >3  in  (15) 
or  (16),  we  find 

2V'2(\'3  +  l)exp|  -  j 
8m(  t)  =  —  1  + - sin 

T 

(17) 


Now  consider  the  limit  of  small  y— 1— >0.  The  limiting 
value  y=  1  corresponds  to  the  value  of  n  =  I  in  (8),  for  which 
the  series  (6)  reduces  to  cos  t.  Expanding  the  coefficients  of 
the  series  (6)  into  a  power  series  near  y=  1 ,  we  obtain  the 
next  term  in  the  asymptotic  expansion  of  8m  in  the  limit  y 
- 1  — >  0,  finite  t  : 


8m(T)  =  cos  t- 


T. 


(19) 


This  formula  describes  oscillations  whose  amplitude 
grows  linearly  with  time.  Matching  (19)  to  the  asymptotic 
expression  (14)  valid  for  large  t,  we  estimate  the  time  when 
8m(r)  reaches  its  peak  as 

Tpeak  —  27/4(y-  1)_1 .  (20) 

Substituting  this  into  (19),  we  obtain  a  rough  estimate  for  the 
peak  value  of  mass  modulation  amplitude 

<Heak=25/4(y—  1)“1/2.  (21) 

Figure  2(c)  compares  the  asymptotic  expression  (19)  to 
the  exact  result  (6)  for  a  very  small  y=  1.003.  The  linear 
growth  (19)  is  seen  to  be  a  good  approximation  for  a  finite 
time.  For  some  time,  before  this  oscillatory  growth  reaches 
the  peak  amplitude  and  evolves  into  an  oscillatory  decay,  the 
growth  is  even  faster  than  Eq.  (19)  predicts. 

The  discussion  of  the  physical  mechanisms  driving  the 
power-law  growth  (16)  and  the  linear  growth  (19)  will  fol¬ 
low  in  Sec.  II  C.  The  next  section  contains  a  rigorous  deri¬ 
vation  of  (16)  for  arbitrary  y. 


B.  Power-law  growth  at  arbitrary  y>  3 

We  start  with  the  continuity  equation  written  in  the  linear 
approximation  for  a  small  single-mode  perturbation  imposed 
on  a  rippled  rarefaction  wave: 

~T  +  ~(pvx)  +  ikp[0)vy  =  0,  (22) 

dt  dx 

where  p(0)  stands  for  the  unperturbed  gas  density  in  the  pla¬ 
nar  rarefaction  wave.  With  the  aid  of  (22)  we  find  that  the 
areal  mass  integral  (2)  satisfies  a  differential  equation 

d  r> 

—  8m  =  -  ikj  p^Vydx,  (23) 


Oscillations  (14)  decay  as  1/r,  the  monotonic  term 
(15) — as  t@,  where 

y—  3 

/8=J— ■  (18) 

y-  1 

For  y<  2  we  have  [3<-  \ ,  and  therefore  the  oscillatory  de¬ 
cay  of  8m  at  large  t  (14)  dominates  over  the  monotonic 
decay  (15).  For  2<  y<  3  the  reverse  is  true. 

In  Figs.  2(a)  and  2(b),  we  compare  the  asymptotic  solu¬ 
tion  (13)  with  the  exact  solution  (6)  at  finite  times  for  mod¬ 
erate  and  large  y.  The  agreement  at  large  r  is  good,  as  it 
should  be. 


where  xj  and  x,  are  the  coordinates  of  the  leading  and  the 
trailing  edges  of  the  expansion  wave,  respectively.  We  nor¬ 
malize  8m  in  the  same  way  as  in  (6),  define  a  dimensionless 
lateral  velocity  8vy  by  vy=-ikSx0asX  8vx,  and  introduce  the 
normalized  speed  of  sound  in  the  unperturbed  expansion 
wave  as  A=a/as.  At  the  leading  edge  of  the  expansion  wave 
A  =  1 ,  at  its  trailing  edge  A  =  0.  The  unperturbed  density  in  the 
expansion  wave  is  expressed  via  A  as  p(°)  =  psA2/(r_1l,  and 
dx=-aastdA,  where  a=(y+  l)/(y-l)  and  ps=ap0  is  the 
postshock  density  [3,4].  Substituting  these  definitions  into 
Eq.  (23),  we  can  present  the  solution  of  this  equation  satis¬ 
fying  the  initial  condition  8m{ 0)=1  as 
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Sm(T )  =  1  -  a2  f  t' dr'  f  A2/(y  ’■* 8uy(A,i,)dA.  (24) 

4+0  J  o 

The  lateral  velocity  profile  Suv(A,t)  is  known  from  the 
linear  perturbation  analysis  of  a  rippled  rarefaction  wave 
[1,4].  For  y>  3  it  is  written  as 

Svv(A,t)  =A^/2  [  cIzSwx(z)Jq(t4>(A,z)),  (25) 

-V 

where  /3  is  given  by  (18),  (f>(A,z)=Aa/2yJn(Al3-z),  n  is  given 
by  (8)  (note  that  for  arbitrary  y>  3,  [3  is  positive,  n  is  nega¬ 
tive,  and  not  necessarily  integer),  and 

2 

8w  |  (z)  =  2  ctjefi EJ_1;  (26) 

j= 0 

2(y-  1)M2  +  y—  3  (y-l)M2-2 

“°  (y-  1)2/?M2  ’  (y-  1)2/?M2  ’ 


y+  1 

a2=  (y-l)2RM2’ 


(27) 


j  =  0,1,2.  (28) 

2  2  p 

Here  M2=(D-U)/as<  1  is  the  downstream  Mach  number  of 
the  incident  shock  wave,  R=D/ (D-U)  >  1  is  the  density 
compression  ratio;  in  the  strong-shock  limit  M2 

— >  \j(y-  l)/(2y),  R— >(y+l)/(y-l)  [see  (3)].  The  function 
8w  \  characterizes  the  initial  profile  of  lateral  velocity  Svy 
produced  when  the  rippled  expansion  wave  was  formed  in 
the  interaction  of  the  planar  incident  shock  wave  with  a 
rippled  free  rear  surface  at  r— >+0  [1,3,4].  Formulas 
(25)-(28)  are  valid  for  a  rippled  rarefaction  wave  produced 
by  a  shock  wave  of  arbitrary  strength  interacting  with  a  ma¬ 
terial  interface  with  arbitrary  density  ratio.  For  y<  3  we  have 
/3<  0,  so  A@>  1,  and  integration  in  (25)  goes  from  1  to  A@. 

One  is  tempted  to  find  the  asymptotic  expression  for 
Sm(r )  at  r— >°°  by  simply  substituting  the  asymptotic  expres¬ 
sion  for  Bessel  function  J0(  T(f>)  at  large  argument  into  (25) 
and  then  performing  the  integrations.  This  approach,  how¬ 
ever,  would  lead  to  incorrect  results:  the  asymptotic  of  a  sum 
is  not  equal  to  the  sum  of  the  asymptotics.  Physically,  we 
have  to  perform  integration  over  standing  sonic  waves  whose 
oscillation  frequencies  are  slightly  different.  When  summed 
up,  the  phase  shifts  between  them  turn  out  to  be  important 
enough  to  drastically  change  the  asymptotic  character  of  the 
result.  Therefore  we  need  a  different  strategy. 

Let  us  substitute  (25)  into  (24)  and  simplify  the  result 
using  the  identity  z,Jo(z)  =  (dl  dz)[_z.J  \(z))  ([12],  Eq.  9.1.30)  to 
integrate  over  r': 


8m(r)  =  1  -  ott2 


f  Aa,2dA  f 
J  0  J  A1 


,  c  ,  Ji (r^(A,z)) 
dz&Wiiz) — — — — 
Ap  t<P(A,z) 


=  1  -  «;£;Q7(t), 


j= o 


(29) 


where 


Q  h)  '■ 


:  f  Aa/2dA  f 

Jo  Ja> 

:-Lf' 

27 71  J, 


,  f._|-/l(7Y/>04,z)) 

dzz  J  - 

ad  t4>(A,z) 


r(-0) 


m 


+100 

d0~  x . 

r(2  +  0)\  4 


2\e 


X  [  dAAan+aeqj(A,0),  j  =  0,1,2.  (30) 

J  o 


In  (30)  we  have  used  integral  representation  ([12],  Eq. 
9.1.26)  for  Ji(z).  The  integrals  over  z  denoted  by  cj}(A ,  6)  in 
(30)  are  expressed  via  hypergeometric  functions  ([13],  Eqs. 
3.197): 


q/A,  9)  ■■ 


■f. 

JAP 


1  dzz**A  1  -  — V  =  A«^)-(1  ^ 


9+  1 


X2F1(9+ei+  1,9+  1;  9+2;  1  -A?). 


(31) 


This  explicit  expression  allows  one  to  perform  the  integra¬ 
tion  over  A  in  (30).  The  result  is 


2/(t)  =  TT 


-J-f 

2/3  J 


, 2a9+ a+2 

+ioo  r  -  e 1  r 

+I  do  \  2/3 


2(3j _ioo  2iri  49+  a+2  +  ^ 

\  2/3  + 

m2^  8 
,  '  4 

x - 7 - 

,  .(  49+  a  +  2 

(9+  1  - +  e, 

\  2(3 


(32) 


We  are  interested  in  the  asymptotic  behavior  of  G7(t")  at 
t  — >  rj~.  The  asymptotic  expressions  sought  for  are  obtained 
by  closing  the  integration  contour  on  the  complex  plane  9  to 
the  left  side  and  calculating  the  residues  contributed  by  the 
poles  at  negative  real  values  of  9.  These  are  zeroes  of  the 
denominator,  9=— 1  and 


°J  =  ~ 


1 

2 


a 

4 


1 

2 


(j  +  a). 


j  =  0,1,2, 


(33) 


and  the  poles  of  the  Gamma  function  in  the  numerator, 

1  1  B 

-  — /,  /  =  0,1,2,  ...  .  (34) 

2  a  a 

After  some  algebra,  we  find  the  asymptotic  expansions  of 
Qj(  t)  in  powers  of  r,  which  for  (29)  translate  into 
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Sm(r)  =  ■ 


~  -  —2j  ajBj 

z  4  j= o 


(  2a0;+  a  +  2 
.  (40,+ a  +  2 

("'+1)rl^-+1 


(-i)T(-e,) 


TjSV 


,_n  ,  /  40;+  a  +  2  \  40,+ a+2 

l~°  Z!(0,+  l)r(  — - +  1  )(  — — - +  e 


2/3 


2/3 


«\6' 


7j(2Z+1)/« 


at  r— >  oo. 


(35) 


Since  for  y>3  we  have  0</3<l  and  «<0,  the  only 
term  on  the  right-hand  side  of  Eq.  (35)  that  grows  at  t— » °o  is 
the  first  j  =  0  term  in  the  square  brackets.  Substituting  (28) 
and  (33)  into  this  term,  and  using  the  duplication  formula 
([12],  Eq.  6.1.18)  to  simplify  the  product  of  two  gamma 
functions  in  the  numerator,  we  transform  it  into 

(y+lV11 

Sm(r)  =  craoIXo')^ - - J  at  t— >  °o.  (36) 

This  is  identical  to  (16).  Validity  of  the  asymptotic  expres¬ 
sion  (16)  for  arbitrary  y>  3  has  thus  been  proved.  Since  the 
y-dependent  coefficient  in  (15)  represents  a  branch  of  the 
same  analytic  function  of  y,  and  it  has  been  checked  in  Sec. 
II  A  that  for  all  y<3  satisfying  (8)  the  corresponding  coef¬ 
ficients  belong  to  the  same  sheet  of  the  Riemann  surface,  we 
conclude  that  (15)  is  also  valid  for  arbitrary  y  between  1  and 
3.  A  similar  argument  applies  to  (14). 

C.  Qualitative  discussion 

The  physical  mechanism  of  unbounded  perturbation 
growth  at  y>  3  is  somewhat  similar  to  that  of  classical  RM 
instability.  In  the  RM  case,  the  fluid  near  a  material  interface 
gets  initial  kick  when  the  incident  planar  shock  wave  inter¬ 
acts  with  it.  During  the  short  time  interval  of  shock-interface 
interaction,  a  finite  surface  vorticity  is  deposited  at  the 
shocked  rippled  interface,  and  this  vorticity  essentially  drives 
the  subsequent  linear  RM  growth. 

Although  the  flow  in  a  rarefaction  wave  propagating  into 
an  unperturbed  gas  is  vorticity-free,  the  shock-free-rippled- 
surface  interaction  produces  a  finite,  constant  in  time  lateral 
fluid  velocity  at  the  trailing  edge  of  the  expansion  wave,  that 
is,  at  the  expansion  front.  Taking  the  limit  A— >0  at  finite  At 
in  (25),  we  find 

Sv V(A  — »  0,  r)  =  -  a0  exp(-  A r\  -  n) .  (37) 

At  the  expansion  front,  where  A  =  0,  we  have  Suy=-a0  at  all 
time.  Of  course,  at  the  trailing  edge  the  density  /i,0>  vanishes. 
However,  since  p-0,=psA2/-7_l7,  and  the  power  2/(y-l)  for 
y>3  is  less  than  unity,  normal  density  gradient  dp^ / dx  di¬ 
verges  at  the  trailing  edge,  indicating  accumulation  of  mass 
in  its  vicinity.  On  the  other  hand,  both  the  unperturbed  pres¬ 
sure  =psA2j,(-7~l]  and  its  normal  gradient  dp^0)/dx  vanish 
at  the  trailing  edge,  which  means  that  mass  could  be  continu¬ 


ously  shifted  laterally  near  the  trailing  edge  without  building 
up  a  substantial  counter-pressure  that  could  give  rise  to  sonic 
oscillations  driving  the  mass  back. 

Figure  3  shows  the  profiles  of  normalized  lateral  velocity 
<%;,,(  77,7)  (a)  and  density  modulation  Sp  [the  latter  is  defined 
by  p= p*01  +  kSx0  X  Sp(y,  r)]  (b)  amplitudes  versus  the  dimen¬ 
sionless  self-similar  coordinate  rj  for  y=5.  Here  r;=(x 
- Ut)/ast=[2-(y+  l)A]/(y-  1)  [3,4];  the  leading  and  the 
trailing  edges  of  the  expansion  wave  correspond  to  77= -1 
and  77=  2/ (y-  1),  respectively.  At  the  trailing  edge,  the  lateral 
velocity  maintains  its  constant  maximum  value  -a0.  Some 
gas  particles,  which  have  initially  received  lateral  velocities 
Suv  close  to  -a0,  continue  their  ballistic  motion  in  the  same 
direction.  If  we  could  assume  that  the  same  set  of  fluid  par¬ 
ticles  has  lateral  velocities,  say,  between  -a0l2  and  -a0  at 
all  times,  then  the  contribution  of  these  particles  to  dm  would 
be  estimated  by  taking  Svy~-a0  and  integrating  over  A  in 
(24)  from  0  to  Amax=constX  T+y-\V(y+7)  (normalized  speed 
of  sound  A  (7)  in  a  given  fluid  particle  in  a  centered  rarefac¬ 
tion  wave  varies  with  time  as  T-Hy-D/fy+i),  see  p],  gq  (17)). 


FIG.  3.  Profiles  of  normalized  lateral  velocity  Svy  (a)  and  den¬ 
sity  modulation  Sp  (b)  amplitudes  versus  self-similar  coordinate  77 
=  (x-Ut)/ast  for  y=5  at  t=  5,  10,  and  15. 
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FIG.  4.  Same  as  in  Fig.  3  for  y=|. 

Not  surprisingly,  this  assumption  translates  into  a  linear 
growth  of  8m  with  time — a  constant  fraction  of  mass  moves 
laterally  with  a  constant  velocity,  resulting  in  a  linear  growth 
of  mass  modulation  amplitude.  This  is  basically  how  the 
classical  RM  instability  develops. 

The  difference  between  the  instability  of  an  expansion 
wave  and  the  classical  RM  instability  is  that  in  the  latter  case 
the  ballistic  fluid  flow  along  the  material  interface  is  vortical 
and,  therefore,  it  is  not  accompanied  with  a  build-up  of  a 
counter-pressure  slowing  down  the  lateral  motion.  On  the 
other  hand,  the  flow  in  the  expansion  wave  is  isoentropic  and 
vorticity-free,  hence  all  our  perturbations  are  sonic  waves. 
Though  the  gas  pressure  near  the  trailing  edge  is  low,  it  is 
still  gradually  slowing  down  the  lateral  motion  of  gas  par¬ 
ticles.  Equation  (37)  shows  that  the  values  of  A  and  t  corre¬ 
sponding  to  Svy  equal  to  some  sizeable  fraction  of  -a0  are 
related  by  A  r=  const.  It  means  that,  as  the  lateral  fluid  motion 
slows  down,  the  number  of  particles  whose  velocities  ex¬ 
ceed,  say,  -a0/2,  decreases  with  time  as  T~2liy+l\  Substitut¬ 
ing  (37)  into  (24)  and  calculating  the  integral  in  the  limit  A 
— * 0  at  finite  At,  we  readily  reproduce  Eq.  (16).  The  density 
profiles  for  this  case  shown  in  Fig.  3(b)  demonstrate  mass 
accumulation  near  the  trailing  edge  of  the  expansion  wave. 

For  comparison.  Fig.  4  shows  the  profiles  of  lateral  veloc¬ 
ity  and  density  in  a  rippled  rarefaction  wave  for  a  stable  case 
of  7=|.  The  lateral  velocity  profiles  of  Fig.  4(a)  are  not  that 
different  from  those  shown  in  Fig.  3(a).  The  gas  mass  near 
the  trailing  edge  is  still  driven  in  the  lateral  direction  at  al¬ 
most  constant  velocity,  but  now  there  is  little  mass  there  to 
redistribute,  and  the  fraction  of  particles  whose  lateral  veloc¬ 
ity  is  close  to  -a0  decreases  so  fast  that  the  total  contribution 
of  this  mechanism  to  8m  decays  with  time.  The  density 
modulation  amplitude  oscillates,  showing  no  accumulation. 

Figure  5  shows  the  profiles  of  lateral  velocity  and  density 
for  y  close  to  unity,  y=  1 . 1 .  Here,  p^  rapidly  decays  with 
increased  distance  from  the  leading  edge  of  the  expansion 


FIG.  5.  Same  as  in  Figs.  3  and  4  for  y=  1.1,  with  profiles  (a)  and 
(b)  zoomed  near  the  head  of  the  expansion  wave,  and  the  profile  of 
8vy  (c)  also  shown  in  full. 

wave.  Therefore,  only  the  gas  layer  near  the  leading  edge 
contributes  to  8m.  Just  as  above,  the  highest  value  of  8uy  is 
near  the  trailing  edge,  but  since  there  is  virtually  no  mass 
there,  ballistic  mass  motion  near  the  trailing  edge  started  in 
the  shock-free-rippled-surface  interaction  adds  no  contribu¬ 
tion  to  8m.  The  areal  mass  is  redistributed  by  the  oscillating 
lateral  velocity  generated  in  the  vicinity  of  the  rippled  lead¬ 
ing  edge  of  the  expansion  wave  as  a  result  of  lateral  pressure 
gradient. 

The  physical  mechanism  of  perturbation  growth  here  is 
essentially  the  same  as  in  a  pendulum  oscillating  in  a  de¬ 
creasing  gravitational  field.  When  the  pendulum  swings  from 
the  equilibrium  position,  the  restoring  force  that  decelerates 
it  is  less  than  the  force  that  accelerated  it  when  it  approached 
the  equilibrium,  hence  in  each  oscillation  it  deviates  farther 
from  the  equilibrium  [14,15].  Similarly,  here  the  pressure 
gradient  drives  the  mass  laterally  until  it  overshoots  the  equi¬ 
librium,  and  a  comparable  counter-pressure  builds  up.  When 
the  gas  mass  starts  moving  in  the  opposite  direction,  the 
unperturbed  pressure  there  is  already  less,  and  it  takes  a 
larger  displacement  to  slow  the  mass  down.  Eventually,  ev¬ 
ery  fluid  particle  approaches  the  trailing  edge  and  the  sonic 
oscillations  start  to  decay  according  to  (14). 

The  same  mechanism  drives  the  so-called  Vishniac  insta¬ 
bility  of  a  blast  wave  in  a  gas  with  low  y  [16].  Low  values  of 
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ideal  gas  y  are  conventionally  used  to  model  dynamics  of 
strongly  radiating  gases.  The  accuracy  of  such  approximation 
is  beyond  the  scope  of  the  present  paper  (see  [17]  and  refer¬ 
ences  therein).  A  blast  wave  consists  of  a  shock  wave  and  a 
rarefaction  wave  immediately  following  it  [18].  A  planar 
shock  wave  in  an  ideal  gas  is  stable  for  any  y  [19].  When 
rippled  shock  and  rarefaction  waves  are  coupled,  as  in  the 
case  of  a  rippled  blast  wave,  the  shock  wave  feeds  to  the 
leading  edge  of  the  rarefaction  wave  a  nonzero  value  of  lat¬ 
eral  velocity  Svy,  in  contrast  with  the  case  of  a  stand-alone 
rippled  rarefaction  wave,  where  8vy  =  0  at  its  leading  edge, 
cf.  Figs.  3(a),  4(a),  and  5(a).  Since  the  rippled  shock  wave 
acts  both  as  a  source  of  perturbations  and  as  a  stabilizer  for 
the  rippled  rarefaction  wave,  all  Vishniac  perturbation 
growth  is  entirely  due  to  the  instability  of  a  rippled  rarefac¬ 
tion  wave  at  low  y  discussed  above.  If  y  is  not  small  enough, 
y>  1.2  for  planar  geometry  [16],  then  there  is  no  instability. 
Even  if  y  is  sufficiently  small,  say,  y=  1 . 1 ,  the  growth  in  a 
planar  geometry  occurs  only  during  a  finite  time  interval  (cf. 
Fig.  5(a)  of  [16]),  same  as  for  a  rippled  expansion  wave  here 
[see  Fig.  1(c)].  A  similar  behavior  of  8m — its  oscillatory 
growth  during  a  finite  time  interval,  which  is  longer  for  a 
lower  y — had  been  established  in  [15]  for  a  shock- 
rarefaction  flow  produced  by  an  impulsive  loading  rather 
than  a  blast. 

Finally,  let  us  discuss  how  the  present  results  are  modified 
for  a  finite-strength  rarefaction  wave,  where  the  normalized 
speed  of  sound  decreases  from  A  =  1  not  to  A  =  0  but  to  some 
finite  minimum  value  of  A  =  M  |  >  0.  This  corresponds  to  the 
unperturbed  shock  wave  incident  upon  a  material  interface 
with  large  density  difference  piight/ Pheavy^  1  rather  than  upon 
a  free  rear/inner  surface.  Then  there  is  no  unlimited  growth 
of  8m  because  all  the  impulsively  accelerated  fluid  particles 
that  would  have  driven  the  growth  in  an  expansion  wave 
leave  the  finite-strength  rarefaction  wave  within  a  finite  time 
interval  r  of  order  of  M[(y+l,/(y_l'1  where 


M,  « 


2ypnBht  y-lV2y 

(  y  —  1 )  Pheavy  ■ 


(38) 


A  finite  counter-pressure  eventually  leads  to  a  phase  reversal 
of  8m,  so  that  at  late  time  the  mass  modulation  amplitude  in 
a  finite  strength  rippled  rarefaction  wave  oscillates  around 
zero  average  value  for  any  y  [see  Fig.  6(a)],  The  growth 
phase  at  y>  3  lasts  only  a  finite  time  and  is  only  noticeable 
for  low  values  of  Mj  that  correspond  to  very  low  density 
ratios  piight/ Pheavy  For  example,  M  [  =0.0 1  in  a  gas  with  y 
=4  corresponds  to  a  strong  shock  interaction  with  a  material 
interface  where  the  density  ratio  is  Plight/ Pheavy =  1-7  X  10-6. 

In  contrast  with  this,  in  a  gas  with  a  low  y  the  behavior  of 
8m  in  a  finite-strength  rarefaction  wave  does  not  differ  from 
the  case  of  an  expansion  wave  [see  the  example  shown  for 
7=1.1  in  Fig.  6(b)].  Here,  the  curves  corresponding  to  M[ 
between  0  and  0.6  are  on  top  of  each  other,  and  the  effect  of 
finite  Mi  is  only  observable  starting  from  M[  =  0.7.  This  is 
because  the  vicinity  of  the  trailing  edge  plays  no  role  here; 
the  mass  redistribution  occurs  predominantly  near  the  lead¬ 
ing  edge,  where  our  exact  solutions  do  not  depend  on  at 
all.  This  is  why  our  qualitative  arguments  are  applicable  to 


FIG.  6.  Time  histories  of  areal  mass  variation  in  a  rippled  rar¬ 
efaction  wave  of  finite  strength,  with  rarefaction  Mach  number 
varying  from  0  (expansion  wave  limit)  to  0.1  for  y=4  (a)  and  to  0.9 
for  7=  1.1  (b). 

the  case  of  a  planar  blast  wave  discussed  above,  where  the 
gas  behind  the  shock  wave  decompresses  to  some  finite, 
time-dependent  minimal  pressure. 

III.  NUMERICAL  SIMULATIONS 

The  theoretical  results  discussed  in  Sec.  II  are  exact,  and 
therefore  do  not  require  numerical  validation.  Rather,  they 
could  be  helpful  for  testing  the  accuracy  of  our  codes  in 
modeling  evolution  of  small  perturbations  in  feedout-related 
problems.  Reproducing  these  exact  theoretical  results  in  a 
numerical  simulation  is  a  challenging  task,  given  the  impor¬ 
tant  roles  of  compressibility,  high  gradients,  and  large  den¬ 
sity  differences  (in  the  range  of  1 04  —  1 06  for  relevant  cases). 
Another  important  factor  is  the  smallness  of  the  initial  am¬ 
plitude  that  ensures  linearity.  If  we  take  the  initial  normalized 
amplitude  kSx0  too  large,  the  nonlinear  effects  will  be  notice¬ 
able  from  the  start.  On  the  other  hand,  small  initial  ampli¬ 
tudes  are  notoriously  hard  to  resolve,  and  there  is  always  a 
risk  either  to  contaminate  the  simulation  results  involving 
such  amplitudes  with  a  numerical  noise  or  to  damp  them 
heavily  with  a  numerical  dissipation.  In  real-life  situations, 
we  have  no  choice  but  to  trust  the  codes;  exact  solutions  help 
us  to  test  how  good  they  are. 

Our  simulations  were  performed  in  two  dimensions  using 
the  FAST2D  hydrocode  developed  at  the  Naval  Research 
Laboratory  [20]  (more  details  and  further  references  are 
given  in  [14])  and  the  Sandbox  code  [21],  which  had  been 
specifically  designed  to  accurately  describe  evolution  of  very 
small  hydrodynamic  perturbations.  The  initial  conditions  cor¬ 
respond  to  the  instant  when  a  planar  shock  wave  reaches  the 
valleys  of  a  rippled  interface  separating  heavy  and  light 
gases.  The  shock  speed  is  D=  107  cm  s"1;  the  ripple  wave¬ 
length  is  \=15  /rm.  Numerical  resolution  in  FAST2D  has 


046306-9 


VELIKOVICH  et  al. 


PHYSICAL  REVIEW  E  72,  046306  (2005) 


FIG.  7.  Simulated  (thick  solid,  dashed,  and  dotted  lines)  and 
exact  analytical  (thin  solid  lines)  time  histories  of  areal  mass  varia¬ 
tion  in  a  rippled  expansion  wave  for  y=  |  and  kSx0= 0.1  (a);  y=4 
and  kSx0  varied  from  2X  10-6  (Sandbox)  to  0.1  to  0.8  (b);  and  y 
=  1.05  and  kSx0= 2X  10“6  (c). 

been  chosen  slightly  higher  than  in  our  modeling  of  the  ex¬ 
periments  on  feedout  [10]  and  is  arguably  adequate  to  the 
practical  state  of  the  art  for  simulations  like  that  performed 
in  ICF-related  studies.  The  version  of  Sandbox  used  here  was 
a  linearized,  fully  spectral  code  in  the  y  direction  and  finite- 
difference  Godunov  hydrodynamic  code  in  the  x  direction. 

We  start  with  y=|,  initial  ripple  amplitude  dx0 
=  0.25  yttm,  and  density  ratio  Plight/ Pheavy =  10~4,  which  corre¬ 
sponds  to  M]=0.24.  Figure  7(a)  demonstrates  a  good  agree¬ 
ment  between  the  theory  and  FAST2D  simulation,  particularly 
in  the  timing  of  minima  and  maxima.  The  remaining  small 
discrepancy  is  not  explained  by  a  finite  value  of  M[  (cf.  Fig. 
6)  because  for  y=|  the  Mj=0  and  M ,  =0.24  theoretical 
curves  are  much  closer  to  each  other  than  theory  and  simu¬ 
lation  curves  in  Fig.  7,  and  small  discrepancy  between  them 
is  seen  both  near  the  minima  and  the  maxima.  To  some  ex¬ 
tent,  the  deviation  seen  in  Fig.  7(a)  may  be  due  to  nonlinear¬ 
ity  (&£>x0  =  0.1),  but  most  likely  it  is  caused  by  numerical 
diffusion. 

A  similar  simulation  has  been  done  for  y=4  and  initial 
ripple  amplitude  <5x0  varied  from  5  X  10  6  (Sandbox)  to 
0.25  to  2  p,m  (fast2d)  [Fig.  7(b)].  The  same  density  ratio 
Plight/ Pheavy =  10' now  corresponds  to  M;  =  0.046.  With  the 
initial  amplitude  dx(}  decreasing,  the  simulated  time  history 


Sm(r)  approaches  the  theoretical  curve,  as  it  should.  Com¬ 
paring  Fig.  7(b)  to  Fig.  6,  we  see  that  the  simulated  curves 
for  &r0=5  X  10-6  p.m  and  0.25  pin  are  actually  closer  to  the 
theoretical  curve  for  M  j  =  0  than  the  actual  theoretical  curve 
for  =  0.046.  Of  course,  the  numerical  dm  includes  contri¬ 
butions  from  the  flow  downstream  from  the  trailing  edge  of 
the  rarefaction  wave,  which,  in  addition  to  numerical  diffu¬ 
sion,  can  slow  down  the  decay  of  dm  with  time  in  compari¬ 
son  to  that  shown  in  Fig.  6. 

Figure  7(c)  presents  the  results  of  the  Sandbox  simulation 
for  y=  1.05,  a  lower  density  ratio  p|lghl/ pheavy=  10-5  that  now 
corresponds  to  Mj  =0.792,  and  initial  ripple  amplitude  dx0 
=  5  X  1(T6  pm.  The  simulation  results  are  in  a  good  agree¬ 
ment  with  the  theoretical  prediction  of  linear  growth,  as  they 
should  be. 


IV.  CONCLUSIONS 

We  have  studied  two  kinds  of  instability  of  a  centered 
expansion  wave.  One  of  them  takes  place  for  large  y>3  and 
manifests  itself  in  an  unlimited  monotonic  power-law  growth 
of  mass  variation  amplitude  with  time  °c t&,  where  [3=(y 
-3)/(y-l).  Of  course,  such  large  values  of  y>3,  while  not 
prohibited  thermodynamically,  do  not  approximate  the  equa¬ 
tion  of  state  of  any  gas  or  other  medium  at  expansion.  Nev¬ 
ertheless,  the  study  of  this  instability,  in  our  opinion,  adds 
some  insight  to  the  general  theory  of  hydrodynamic  instabili¬ 
ties.  The  power  index  of  this  new  instability  is  determined  by 
the  equation  of  state,  just  as  for  Vishniac’s  instability  of  a 
blast  wave,  and,  in  contrast  with  the  latter  case,  is  expressed 
in  a  closed  analytical  form.  The  physical  mechanism  driving 
this  instability,  however,  is  entirely  different  from  Vishniac’s 
instability.  Rather,  it  is  essentially  the  same  as  in  classical 
Richtmyer-Meshkov  instability.  The  instability  is  driven  by  a 
ballistic  motion  of  the  fluid  in  the  lateral  direction,  which  is 
launched  by  the  initial  interaction  of  a  smooth  incident  shock 
with  a  rippled  rear  surface.  We  can  state  that  this  new  insta¬ 
bility  is  in  a  sense  an  expansion-wave  counterpart  of  the 
shock-induced  classical  RM  instability.  In  the  RM  case,  the 
fluid  initially  put  into  the  lateral  motion  continues  to  slide 
along  the  material  interface  without  creating  a  counter¬ 
pressure  that  could  slow  down  its  flow,  which  results  in  a 
linear  growth  of  mass  variation  amplitude  dm  with  time.  In 
an  expansion  wave,  the  lateral  flow  of  a  fluid  is  a  sonic  wave 
in  a  gas  whose  speed  of  sound  in  the  vicinity  of  the  expan¬ 
sion  front  approaches  but  not  equals  zero.  The  counter¬ 
pressure  developing  in  such  a  sonic  wave  gradually  slows 
down  the  fluid  particles,  removing  them  from  the  initially 
launched  lateral  flow.  Our  condition  y>  3  is  required  to 
make  the  contribution  of  these  particles  to  dm  grow  faster 
than  the  rate  at  which  they  are  removed.  If  this  condition  is 
satisfied,  then  the  power  index  f3  is  still  less  than  unity,  the 
power  of  time  characteristic  of  the  classical  RM  instability, 
as  it  should  be. 

We  have  also  studied  another  kind  of  instability  which 
develops  for  y  close  to  unity.  It  is  sometimes  argued  that 
such  values  of  y  represent  reasonable,  if  crude,  approxima¬ 
tion  of  equation  of  state  for  strongly  radiating  gases.  This 
instability  develops  not  only  in  expansion  waves  but  in  rar- 
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efaction  waves  of  arbitrary  strength.  It  manifests  itself  in  an 
oscillatory  growth  of  Sin ,  the  oscillations  amplitude  increas¬ 
ing  linearly  with  time.  The  growth  continues  for  a  finite  time 
interval,  which  is  longer  for  y  closer  to  unity,  until  dm 
reaches  a  peak  value,  which  is  higher  for  y  closer  to  unity, 
and  then  an  oscillatory  decay  follows.  Similar  behavior  had 
been  demonstrated  for  Vishniac’s  instability  of  a  blast  wave 
in  planar  geometry  [16].  In  this  case,  apparent  similarity  with 
Vishniac’s  instability  is  based  on  a  physical  analogy.  In  our 
opinion,  Vishniac’s  perturbation  growth  in  a  blast  wave  is 
driven  by  exactly  the  same  mechanism  as  described  here,  the 
instability  developing  in  the  rarefaction-wave  component  of 
a  rippled  blast  wave.  Its  shock-wave  component  serves  both 


as  a  stabilizer  and  as  a  source  of  initial  perturbations  for  the 
rippled  rarefaction  wave  that  follows  it  in  a  blast  wave,  but  is 
not  unstable  by  itself  for  any  y. 
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